// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtAMRFileFormat.C                             //
// ************************************************************************* //
#include <avtAMRFileFormat.h>

#include <vtkCellType.h>
#include <vtkFieldData.h>
#include <vtkFloatArray.h>
#include <vtkIntArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>
#include <avtIntervalTree.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtStructuredDomainNesting.h>
#include <avtVariableCache.h>

#include <DebugStream.h>
#include <Expression.h>
#include <InvalidDBTypeException.h>
#include <InvalidVariableException.h>


#include <AMRreaderWithLevels.h>

#include <string>
#include <stdlib.h>

#ifdef ENGINE
#define DO_DOMAIN_NESTING
#define DO_DOMAIN_BOUNDARIES
#endif

#define DEBUG_PRINT

const std::string avtAMRFileFormat::amr_name     ="AMR";
const std::string avtAMRFileFormat::intf_name    ="Interface";
const std::string avtAMRFileFormat::invalid_name ="Invalid";

std::string avtAMRFileFormat::
composeName( const std::string& m, const std::string& v, const char app )
{
    return m+app+v;
}


void avtAMRFileFormat::
decomposeName( const std::string& s, std::string& m, std::string& v )
{
    if( s.compare( 0, amr_name.size(), amr_name )==0 )
        m = amr_name;
    else if( s.compare( 0, intf_name.size(), intf_name )==0 )
        m = intf_name;
    else
        m = invalid_name;

    if( m!=invalid_name )
        v = s.substr( m.size()+1 );
    else
        v = s;
}



int avtAMRFileFormat::
GetCycle()
{
    return GetReader(false)->GetNumCycles();
}


double avtAMRFileFormat::
GetTime()
{
    return GetReader(false)->GetSimuTime();
}



// ****************************************************************************
//  Method: avtAMRFileFormat constructor
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
// Modifications:
//   Brad Whitlock, Wed May 28 16:15:20 PDT 2014
//   Change how the reader is instantiated. Add enableAMR to indicate whether
//   new AMR enhancements will be enabled.
//
// ****************************************************************************

avtAMRFileFormat::avtAMRFileFormat(const char *filename)
    : avtSTMDFileFormat(&filename, 1)
{
    reader_ = NULL;

    // Enable the new features by default unless we're purposefully disabling it.
    enableAMR = getenv("AMRDISABLENEW")==NULL;

    // Check the file only if we're in strict mode. This will prevent many files
    // from being opened upfront.
    if(this->GetStrictMode())
        GetReader(false);
}


avtAMRFileFormat::~avtAMRFileFormat()
{
    debug2 << "closed avtAMRFileFormat\n.";
    if(reader_ != NULL)
    {
        delete reader_;
    }
}

// ****************************************************************************
// Method: avtAMRFileFormat::GetReader
//
// Purpose:
//   Get the reader object.
//
// Note:       We can defer creation of the real object until we access it from
//             calling this method.
//
// Programmer: Brad Whitlock
// Creation:   Wed May 28 16:13:58 PDT 2014
//
// Modifications:
//
// ****************************************************************************

AMRreaderInterface *
avtAMRFileFormat::GetReader(bool heavyweightInit)
{
    if(reader_ == NULL)
    {
        // Create the reader that will read the file.
        if(enableAMR)
        {
            reader_ = new AMRreaderWithLevels;
        }
        else
        {
            reader_ = InstantiateAMRreader();
        }

        // Read the file information.
        int ierr = reader_->getInfo(filenames[0]);
        if( ierr!=0 )
        {
            EXCEPTION1(InvalidDBTypeException,
                       "The file could not be opened");
        }
    }

    // Do more time-consuming stuff only when we need to.
    if(heavyweightInit)
    {
        reader_->BuildMetaData();
    }

    return reader_;
}


// ****************************************************************************
//  Method: avtAMRFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
// ****************************************************************************

void
avtAMRFileFormat::FreeUpResources(void)
{
    debug2 << "Now free resource.\n";
    if(reader_ != NULL)
    {
        //reader_->freedata();
        delete reader_;
        reader_ = NULL;
    }
    debug2 << "Now done free resource.\n";
}


// ****************************************************************************
//  Method: avtAMRFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
//  Modifications:
//    Brad Whitlock, Tue May 27 11:10:58 PDT 2014
//    I changed the AMR mesh into an AMR mesh. I also added domain nesting and
//    domain boundaries so VisIt can ghost AMR levels and remove internal faces
//    within levels.
//
// ****************************************************************************

void
avtAMRFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    // meshes
    int nblks = GetReader()->GetNumberOfBlocks();
    if( nblks>0 )
    {
        // Basic mesh metadata
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = amr_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 3;
        mmd->meshType = enableAMR ? AVT_AMR_MESH : AVT_RECTILINEAR_MESH;
        mmd->numBlocks = nblks;

        if(enableAMR)
        {
            int nLevels = GetReader()->GetNumberOfLevels();
#ifdef DEBUG_PRINT
            debug1 << "Num levels: " << nLevels << endl;
#endif
            // Create ratios for avtStructuredDomainNesting.
#ifdef DO_DOMAIN_NESTING
            avtStructuredDomainNesting *dn = new avtStructuredDomainNesting(nblks, nLevels);
            dn->SetNumDimensions(3);
            bool *levelCSset = new bool[nLevels];
            for(int level = 0; level < nLevels; ++level)
            {
                std::vector<int> ratios;
                int r = (level == 0) ? 1 : 2;
                ratios.push_back(r);
                ratios.push_back(r);
                ratios.push_back(r);
                dn->SetLevelRefinementRatios(level, ratios);
                levelCSset[level] = false;
            }
#endif
#ifdef DO_DOMAIN_BOUNDARIES
            // Create ratios for avtRectilinearDomainBoundaries
            avtStructuredDomainBoundaries *sdb = new avtRectilinearDomainBoundaries(true);
            sdb->SetNumDomains(nblks);
#endif

            // Use information about each block to fill in metadata, domain nesting,
            // and domain boundaries.
            std::vector<int>         groupIds(nblks);
            std::vector<std::string> blockPieceNames(nblks);
            int localPatchNumber[20];
            memset(localPatchNumber, 0, sizeof(int) * 20);
            for(int bid = 0; bid < nblks; ++bid)
            {
                // Get the level and level-appropriate IJK indexing for the block.
                int level = 0, ijk_start[3] = {0,0,0}, ijk_end[3] = {0,0,0};
                GetReader()->GetBlockHierarchicalIndices(bid, &level, ijk_start, ijk_end);

                //
                // Fill in some fields for the metadata.
                // Use the level and local patch number to make up a name for the patch.
                //
                // COMMENT: If this naming scheme introduces too much shifting
                //          from time step to time step, maybe identifying the patch
                //          using a representation of its block key could be better.
                //
                char tmpName[128];
                groupIds[bid] = level;
                snprintf(tmpName, 128, "level%02d,patch%04d", level, localPatchNumber[level]++);
                blockPieceNames[bid] = tmpName;
#ifdef DEBUG_PRINT
                debug1 << "Block " << bid << " level: " << groupIds[bid]
                       << ", name: " << blockPieceNames[bid] << endl;
#endif
                //
                // Fill in avtStructuredDomainNesting.
                //
#ifdef DO_DOMAIN_NESTING
                if(!levelCSset[level])
                {
                    float xs[3], dx[3];
                    GetReader()->GetBlockMesh(bid, xs, dx);
                    std::vector<double> cs;
                    cs.push_back(dx[0]);
                    cs.push_back(dx[1]);
                    cs.push_back(dx[2]);
                    dn->SetLevelCellSizes(level, cs);
                    levelCSset[level] = true;
                }

                std::vector<int> logicalExtents(6);
                logicalExtents[0] = ijk_start[0];
                logicalExtents[1] = ijk_start[1];
                logicalExtents[2] = ijk_start[2];
                logicalExtents[3] = ijk_end[0];
                logicalExtents[4] = ijk_end[1];
                logicalExtents[5] = ijk_end[2];
#ifdef DEBUG_PRINT
                debug1 << "Block " << bid << " logical indices: "
                       << logicalExtents[0] << ", "
                       << logicalExtents[1] << ", "
                       << logicalExtents[2] << ", "
                       << logicalExtents[3] << ", "
                       << logicalExtents[4] << ", "
                       << logicalExtents[5] << endl;
#endif
                // Get the child patches. The patches are not necessarily ordered in the list but
                // all patches (except 0) will have a direct parent.
                std::vector<int> childPatches;
                if(bid == 0)
                {
                    for(int cbid = 1; cbid < nblks; ++cbid)
                    {
                        if(OctKey_NumLevels(GetReader()->GetBlockKey(cbid)) == 2)
                            childPatches.push_back(cbid);
                    }
                }
                else
                {
                    OctKey thisKey(GetReader()->GetBlockKey(bid));
                    for(int cbid = 1; cbid < nblks; ++cbid)
                    {
                        if(cbid != bid && OctKey_HasImmediateParent(GetReader()->GetBlockKey(cbid), thisKey))
                            childPatches.push_back(cbid);
                    }
                }
#ifdef DEBUG_PRINT
                debug1 << "Block " << bid << " child patches: ";
                for(size_t q = 0; q < childPatches.size(); ++q)
                    debug1 << childPatches[q] << ", ";
                debug1 << endl;
#endif
                dn->SetNestingForDomain(bid, level, childPatches, logicalExtents);
#endif

                //
                // Fill in avtRectilinearDomainBoundaries.
                //
#ifdef DO_DOMAIN_BOUNDARIES
                int e[6];
                e[0] = ijk_start[0];
                e[1] = ijk_end[0] + 1;
                e[2] = ijk_start[1];
                e[3] = ijk_end[1] + 1;
                e[4] = ijk_start[2];
                e[5] = ijk_end[2] + 1;
                sdb->SetIndicesForAMRPatch(bid, level, e);
#endif
            }

#ifdef DO_DOMAIN_NESTING
            delete [] levelCSset;
            debug2 << "Caching domain nesting." << endl;
            void_ref_ptr vr = void_ref_ptr(dn, avtStructuredDomainNesting::Destruct);
            cache->CacheVoidRef(amr_name.c_str(), AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                                timestep, -1, vr);
#endif
#ifdef DO_DOMAIN_BOUNDARIES
            debug2 << "Caching domain boundaries." << endl;
            sdb->CalculateBoundaries();
            void_ref_ptr vsdb = void_ref_ptr(sdb,avtStructuredDomainBoundaries::Destruct);
            // HACK: VisIt won't use this information unless "any_mesh" is used.
            cache->CacheVoidRef("any_mesh", AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                                timestep, -1, vsdb);
#endif

            // AMR mesh metadata
            mmd->blockTitle = "patches";
            mmd->blockPieceName = "patch";
            mmd->blockNames = blockPieceNames;

            mmd->numGroups = nLevels;
            mmd->groupTitle = "levels";
            mmd->groupPieceName = "level";
            mmd->groupIds = groupIds;
        }

        // Add the mesh metadata.
        md->Add(mmd);
    }

    int np,ne;
    GetReader()->GetInterfaceSizes( &np, &ne );
    if( np>0 && ne>0 )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = intf_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 2;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = 1;
        md->Add(mmd);
    }

    // AMR variable
    if( nblks>0 )
    {
        // scalar
        avtScalarMetaData *smd;

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "density");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "u velocity");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "v velocity");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "w velocity");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "pressure");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "temperature");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "sound speed");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "internal energy");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "kinetic energy");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( amr_name, "level");
        smd->meshName = amr_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        // vector
        avtVectorMetaData *vmd;
        vmd = new avtVectorMetaData;
        vmd->name = composeName( amr_name, "velocity");
        vmd->meshName = amr_name;
        vmd->centering = AVT_ZONECENT;
        vmd->hasUnits = false;
        // ??? is it necessary
        //vmd->varDim = 3;
        md->Add(vmd);

        //tag
        if( GetReader()->HasTag() )
        {
            smd = new avtScalarMetaData;
            smd->name = composeName( amr_name, "cell tag");
            smd->meshName = amr_name;
            smd->centering = AVT_ZONECENT;
            smd->hasUnits = false;
            md->Add(smd);
        }
    }



    // interface
    if( np>0 && ne>0 )
    {
        // scalar
        avtScalarMetaData *smd;

        smd = new avtScalarMetaData;
        smd->name = composeName( intf_name, "u velocity");
        smd->meshName = intf_name;
        smd->centering = AVT_NODECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( intf_name, "v velocity");
        smd->meshName = intf_name;
        smd->centering = AVT_NODECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( intf_name, "w velocity");
        smd->meshName = intf_name;
        smd->centering = AVT_NODECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( intf_name, "pressure");
        smd->meshName = intf_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        smd = new avtScalarMetaData;
        smd->name = composeName( intf_name, "material");
        smd->meshName = intf_name;
        smd->centering = AVT_ZONECENT;
        smd->hasUnits = false;
        md->Add(smd);

        // vector
        avtVectorMetaData *vmd;
        vmd = new avtVectorMetaData;
        vmd->name = composeName( intf_name, "velocity");
        vmd->meshName = intf_name;
        vmd->centering = AVT_NODECENT;
        vmd->hasUnits = false;
        // ??? is it necessary
        //vmd->varDim = 3;
        md->Add(vmd);
    }
}


// ****************************************************************************
//  Method: avtAMRFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
// ****************************************************************************

vtkDataSet *
avtAMRFileFormat::GetMesh(int domain, const char *meshname)
{
    // amr
    if( amr_name.compare(meshname)==0 )
    {
        int blkdim[3];
        GetReader()->GetBlockDimensions( domain, blkdim );
        // count the # of points
        ++blkdim[0];
        ++blkdim[1];
        ++blkdim[2];

        float xs[3],dx[3];
        int ierr = GetReader()->GetBlockMesh( domain, xs, dx );
        if( ierr!=0 )
        {
            EXCEPTION1( InvalidVariableException,
                        "Failed to get given block mesh.\n" );
        }

        vtkFloatArray *crd[3];
        for( int as=0; as<3; as++ )
        {
            crd[as] = vtkFloatArray::New();
            crd[as]->SetNumberOfTuples( blkdim[as] );
            float* vv = (float*)crd[as]->GetVoidPointer(0);
            for( int i=0; i<blkdim[as]; i++ )
                vv[i] = xs[as] + (float)i*dx[as];
        }

        vtkRectilinearGrid *grd = vtkRectilinearGrid::New();
        grd->SetDimensions( blkdim );
        grd->SetXCoordinates( crd[0] );
        grd->SetYCoordinates( crd[1] );
        grd->SetZCoordinates( crd[2] );

        for( int as=0; as<3; as++ )
            crd[as]->Delete();

#if 0
        // Needed for AMR Stitch Cell support...
        if(enableAMR)
        {
            // Get the level and level-appropriate IJK indexing for the block.
            int level = 0, ijk_start[3] = {0,0,0}, ijk_end[3] = {0,0,0};
            GetReader()->GetBlockHierarchicalIndices(domain, &level, ijk_start, ijk_end);

            vtkIntArray *arr = vtkIntArray::New();
            arr->SetNumberOfTuples(3);
            arr->SetValue(0, ijk_start[0]);
            arr->SetValue(1, ijk_start[1]);
            arr->SetValue(2, ijk_start[2]);
            arr->SetName("base_index");
            grd->GetFieldData()->AddArray(arr);
            arr->Delete();
        }
#endif

        return grd;
    }

    // interface
    if( intf_name.compare(meshname)==0 )
    {
        int np,ne;
        GetReader()->GetInterfaceSizes( &np, &ne );

        if(np <= 0 || ne <= 0)
        {
            debug2 << "InterfaceSizes returned bad values: np=" << np << ", ne=" << ne << endl;
            return NULL;
        }

        vtkPoints* crd = vtkPoints::New();
        crd->SetNumberOfPoints(np);
        float *pts = (float *) crd->GetVoidPointer(0);
        int ierr = GetReader()->GetInterfaceVariable( AMRreaderInterface::i_coor, pts );
        if( ierr!=0 )
        {
            crd->Delete();
            EXCEPTION1( InvalidVariableException,
                        "Failed to read in interface coordinates.\n" );
        }

        int *elmt = new int[ 4*ne ];
        ierr = GetReader()->GetInterfaceVariable( AMRreaderInterface::i_segt, elmt );
        if( ierr!=0 )
        {
            crd->Delete();
            delete [] elmt;
            std::string msg="Failed to read in interface elements of ";
            msg+=meshname;
            EXCEPTION1( InvalidVariableException, msg );
        }

        vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
        ugrid->SetPoints( crd );
        ugrid->Allocate( ne );

        vtkIdType verts[4];
        int* cnt=elmt;
        for( int i=0; i<ne; i++ )
        {
            verts[0] = *cnt++;
            verts[0]-=1;
            verts[1] = *cnt++;
            verts[1]-=1;
            verts[2] = *cnt++;
            verts[2]-=1;
            verts[3] = *cnt++;
            verts[3]-=1;
            ugrid->InsertNextCell( VTK_QUAD, 4, verts );
        }

        delete [] elmt;

        return ugrid;
    }

    std::string msg="Unrecognized mesh name: ";
    msg+=meshname;
    EXCEPTION1( InvalidVariableException, msg );
}


// ****************************************************************************
//  Method: avtAMRFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
// Modifications:
//   Brad Whitlock, Fri Jun  6 11:15:50 PDT 2014
//   Free objects before throwing exceptions.
//
// ****************************************************************************

vtkDataArray *
avtAMRFileFormat::GetVar(int domain, const char *name)
{
    std::string mesh,varname;
    decomposeName( name, mesh, varname );
    debug2 << "geting variable " << varname << " on mesh " << mesh << "...\n";

    // amr
    if( mesh.compare(amr_name)==0 )
    {
        int vid;
        if( varname.compare("density")==0 )
            vid=AMRreaderInterface::v_dens;
        else if( varname.compare("u velocity")==0 )
            vid=AMRreaderInterface::v_uvel;
        else if( varname.compare("v velocity")==0 )
            vid=AMRreaderInterface::v_vvel;
        else if( varname.compare("w velocity")==0 )
            vid=AMRreaderInterface::v_wvel;
        else if( varname.compare("pressure")==0 )
            vid=AMRreaderInterface::v_pres;
        else if( varname.compare("temperature")==0 )
            vid=AMRreaderInterface::v_temp;
        else if( varname.compare("sound speed")==0 )
            vid=AMRreaderInterface::v_sndv;
        else if( varname.compare("internal energy")==0 )
            vid=AMRreaderInterface::v_eint;
        else if( varname.compare( "kinetic energy")==0 )
            vid=AMRreaderInterface::v_eknt;
        else if( varname.compare( "cell tag")==0 )
            vid=AMRreaderInterface::v_tags;
#if 1
        else if( enableAMR && varname.compare( "level")==0 )
        {
            int sz = GetReader()->GetBlockSize( domain );
            vtkFloatArray* var = vtkFloatArray::New();
            var->SetNumberOfTuples( sz );

            OctKey k = GetReader()->GetBlockKey(domain);
            int level = OctKey_NumLevels(k)-1;
            for(int i = 0; i < sz; ++i)
                var->SetTuple1(i, level);
            return var;
        }
#endif
        else
        {
            std::string msg = "Unknown scalar variable ";
            msg+=name;
            msg+=" is specified.";
            EXCEPTION1( InvalidVariableException, msg );
        }

        int sz = GetReader()->GetBlockSize( domain );
        debug5 << "vid=" << vid << "  domain=" << domain << "  sz=" << sz << "\n";
        vtkFloatArray* var = vtkFloatArray::New();
        var->SetNumberOfTuples( sz );
        float *val = (float *)var->GetVoidPointer(0);

        int err = GetReader()->GetBlockVariable( domain, vid, val );
        if( err!=0 )
        {
            var->Delete();
            std::string msg = "Failed to read in scalar variable ";
            msg+=name;
            EXCEPTION1( InvalidVariableException, msg );
        }
        return var;
    }

    // interface
    if( mesh.compare(intf_name)==0 )
    {
        int np,ne;
        GetReader()->GetInterfaceSizes( &np, &ne );

        int vid=0,cmp=0,sz=0;
        if( varname.compare("u velocity")==0 )
        {
            vid=AMRreaderInterface::i_velo;
            cmp=0;
            sz=np;
        }
        else if( varname.compare("v velocity")==0 )
        {
            vid=AMRreaderInterface::i_velo;
            cmp=1;
            sz=np;
        }
        else if( varname.compare("w velocity")==0 )
        {
            vid=AMRreaderInterface::i_velo;
            cmp=2;
            sz=np;
        }
        else if( varname.compare("pressure")==0 )
        {
            vid=AMRreaderInterface::i_pres;
            sz=ne;
        }
        else if( varname.compare("material")==0 )
        {
            vid=AMRreaderInterface::i_matl;
            sz=ne;
        }
        else
        {
            std::string msg = "Unknown scalar variable ";
            msg+=name;
            msg+=" is specified.";
            EXCEPTION1( InvalidVariableException, msg );
        }

        vtkFloatArray* var = vtkFloatArray::New();
        var->SetNumberOfTuples( sz );
        float *val = (float *)var->GetVoidPointer(0);

        if( vid==AMRreaderInterface::i_velo )
        {
            float* buf = new float[sz*3];
            int err = GetReader()->GetInterfaceVariable( vid, buf );
            if( err!=0 )
            {
                delete [] buf;
                var->Delete();
                std::string msg = "Failed to read in scalar variable ";
                msg+=name;
                EXCEPTION1( InvalidVariableException, msg );
            }

            for( int i=0; i<np; i++ )
                val[i] = buf[3*i+cmp];

            delete [] buf;
        }
        else
        {
            int err = GetReader()->GetInterfaceVariable( vid, val );
            if( err!=0 )
            {
                var->Delete();
                std::string msg = "Failed to read in scalar variable ";
                msg+=name;
                EXCEPTION1( InvalidVariableException, msg );
            }
        }
        return var;
    }

    std::string msg="Unrecongnized mesh for scalar variable ";
    msg+=name;
    EXCEPTION1( InvalidVariableException, msg );
}


// ****************************************************************************
//  Method: avtAMRFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Apr 26 10:26:32 PDT 2010
//
// Modifications:
//   Brad Whitlock, Fri Jun  6 11:15:50 PDT 2014
//   Free objects before throwing exceptions.
//
// ****************************************************************************

vtkDataArray *
avtAMRFileFormat::GetVectorVar(int domain, const char *name)
{
    std::string mesh,varname;
    decomposeName( name, mesh, varname );
    debug2 << "geting variable " << varname << " on mesh " << mesh << "...\n";

    if( mesh.compare(amr_name)==0 )
    {
        int vid;
        if( varname.compare("velocity")==0 )
            vid=AMRreaderInterface::v_velo;
        else
        {
            std::string msg = "Unknown vector variable ";
            msg+=name;
            msg+=" is specified.";
            EXCEPTION1( InvalidVariableException, msg );
        }

        int sz = GetReader()->GetBlockSize( domain );
        vtkFloatArray* var = vtkFloatArray::New();
        var->SetNumberOfComponents(3);
        var->SetNumberOfTuples( sz );
        float *val = (float *)var->GetVoidPointer(0);

        int err = GetReader()->GetBlockVariable( domain, vid, val );
        if( err!=0 )
        {
            var->Delete();
            std::string msg="Failed to read in vector variable ";
            msg+=name;
            EXCEPTION1( InvalidVariableException, msg );
        }
        return var;
    }


    if( mesh.compare(intf_name)==0 )
    {
        int np,ne;
        GetReader()->GetInterfaceSizes( &np, &ne );

        int vid;
        if( varname.compare("velocity")==0 )
            vid=AMRreaderInterface::i_velo;
        else
        {
            std::string msg = "Unknown vector variable ";
            msg+=name;
            msg+=" is specified.";
            EXCEPTION1( InvalidVariableException, msg );
        }

        vtkFloatArray* var = vtkFloatArray::New();
        var->SetNumberOfComponents(3);
        var->SetNumberOfTuples( np );
        float *val = (float *)var->GetVoidPointer(0);

        int err = GetReader()->GetInterfaceVariable( vid, val );
        if( err!=0 )
        {
            var->Delete();
            std::string msg="Failed to read in vector variable ";
            msg+=name;
            EXCEPTION1( InvalidVariableException, msg );
        }
        return var;
    }

    std::string msg="Unrecognized mesh for vector variable ";
    msg+=name;
    EXCEPTION1( InvalidVariableException, msg );
}

// ****************************************************************************
// Method: avtAMRFileFormat::GetAuxiliaryData
//
// Purpose:
//   Return spatial metadata.
//
// Arguments:
//
//
// Returns:
//
// Note:
//
// Programmer: Brad Whitlock
// Creation:   Fri May 23 10:53:20 PDT 2014
//
// Modifications:
//
// ****************************************************************************

void *
avtAMRFileFormat::GetAuxiliaryData(const char *var, int domain,
                                   const char *type, void *args, DestructorFunction &df)
{
    void *retval = 0;

    if(enableAMR && amr_name.compare(var) == 0)
    {
        if(strcmp(type, AUXILIARY_DATA_SPATIAL_EXTENTS) == 0)
        {
            avtIntervalTree *itree = new avtIntervalTree(GetReader()->GetNumberOfBlocks(), 3);
            for(int dom = 0; dom < GetReader()->GetNumberOfBlocks(); ++dom)
            {
                int blkdim[3]= {0,0,0};
                GetReader()->GetBlockDimensions(dom, blkdim);

                float xs[3]= {0.f,0.f,0.f}, dx[3]= {0.f,0.f,0.f};
                GetReader()->GetBlockMesh(dom, xs, dx);

                double extents[6];
                extents[0] = xs[0];
                extents[1] = xs[0] + blkdim[0] * dx[0];
                extents[2] = xs[1];
                extents[3] = xs[1] + blkdim[1] * dx[1];
                extents[4] = xs[2];
                extents[5] = xs[2] + blkdim[2] * dx[2];

                itree->AddElement(dom, extents);
            }
            itree->Calculate(true);
            retval = (void *)itree;
            df = avtIntervalTree::Destruct;
        }
    }

    return retval;

}
